Environmental Data Preperation

Coral Growth Dataset Structure

The coral growth data for this analysis variables related to the growth of individual colonies collected from the visual analysis of x-rays of 1 inch wide slices cut down the central growth plane of individual colonies from three Persian Gulf sites and two sites off the coast of Oman. The growth variables are measured in yearly increments, so environmental data must also be summarized in yearly increments.

Lets look at the structure of the growth data first:

onyr <- read.csv("rawdata/maxvalues-year2000onwards.csv")
head(onyr)

##   year   site colony track   sample  calc extension  density hdbands partmort
## 1 2000 AlAqah  ala12    t1 ala12_t1 0.797    0.6559 1.534612       0        0
## 2 2000 AlAqah  ala12    t2 ala12_t2 1.007    0.5701 1.397308       0        0
## 3 2000  Dibba  dib01    t1 dib01_t1 0.956    0.7467 1.280707       0        0
## 4 2000  Dibba  dib01    t2 dib01_t2 0.526    0.3784 1.390043       0        0
## 5 2000  Dibba  dib04    t1 dib04_t1 0.940    0.8173 1.150707       0        0
## 6 2000  Dibba  dib04    t2 dib04_t2 1.335    1.0998 1.213488       0        0
##   height dia
## 1   13.5  23
## 2   13.5  23
## 3     NA  NA
## 4     NA  NA
## 5   14.0  28
## 6   14.0  28

summary(onyr)

##       year          site              colony             track          
##  Min.   :2000   Length:768         Length:768         Length:768        
##  1st Qu.:2005   Class :character   Class :character   Class :character  
##  Median :2008   Mode  :character   Mode  :character   Mode  :character  
##  Mean   :2008                                                           
##  3rd Qu.:2011                                                           
##  Max.   :2013                                                           
##                                                                         
##     sample               calc          extension         density      
##  Length:768         Min.   :0.1650   Min.   :0.1665   Min.   :0.8693  
##  Class :character   1st Qu.:0.5965   1st Qu.:0.4440   1st Qu.:1.1948  
##  Mode  :character   Median :0.7695   Median :0.5750   Median :1.3423  
##                     Mean   :0.8211   Mean   :0.6044   Mean   :1.3722  
##                     3rd Qu.:1.0170   3rd Qu.:0.7232   3rd Qu.:1.5187  
##                     Max.   :1.9930   Max.   :2.9767   Max.   :2.3394  
##                     NA's   :4        NA's   :9        NA's   :7       
##     hdbands           partmort            height            dia       
##  Min.   :0.00000   Min.   :0.000000   Min.   : 4.000   Min.   :11.00  
##  1st Qu.:0.00000   1st Qu.:0.000000   1st Qu.: 7.000   1st Qu.:15.00  
##  Median :0.00000   Median :0.000000   Median : 8.500   Median :16.00  
##  Mean   :0.08333   Mean   :0.003906   Mean   : 8.814   Mean   :17.94  
##  3rd Qu.:0.00000   3rd Qu.:0.000000   3rd Qu.:11.000   3rd Qu.:22.00  
##  Max.   :1.00000   Max.   :1.000000   Max.   :14.000   Max.   :28.00  
##                                       NA's   :50       NA's   :50

# number of unique colonies
length(unique(onyr$colony))

## [1] 38

# number of unique sites = 5
unique(onyr$site) 

## [1] "AlAqah"     "Dibba"      "RasGhanada" "Saadiyat"   "Delma"

# plot of individual colony diameter Vs height
xyplot(height ~ dia | site, data = onyr,
       #group = year, #t="l",
       grid = TRUE, auto.key=list(space="top", columns=5, 
                                  title="Site", cex.title=1,
                                  lines=TRUE, points=FALSE))

So we have 768 measurements from 38 colonies from five sites. The coral growth data record growth from 2000-2013.

Environmental Variable Summaries

Lets look at the structure of the environmental variables, which are data products derived from satellite observations. They consist of sea surface temperature metrics and current speed and direction metrics from each of the sites.

Here is the Sea Surface Temperature (SST) data for each site. We can see that the data product summarizes the SST daily and the dataset starts in 1985:

temp <- read.csv("rawdata/ctemp_sst_mod.csv")
str(temp)

## 'data.frame':    12370 obs. of  11 variables:
##  $ date      : chr  "1985JAN01" "1985JAN02" "1985JAN03" "1985JAN04" ...
##  $ date1     : chr  "1/1/85" "2/1/85" "3/1/85" "4/1/85" ...
##  $ yearday   : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Year      : int  1985 1985 1985 1985 1985 1985 1985 1985 1985 1985 ...
##  $ Delma     : num  22.7 22.6 22.5 22.2 22.1 ...
##  $ RasGhanada: num  23.3 23.1 23 22.8 22.8 ...
##  $ Saadiyat  : num  23.3 23.1 23 22.7 22.7 ...
##  $ AlAqah    : num  24 23.9 23.9 23.9 23.9 ...
##  $ Dibba     : num  24 23.9 23.9 23.8 23.8 ...
##  $ month     : chr  "JAN" "JAN" "JAN" "JAN" ...
##  $ season    : chr  "winter" "winter" "winter" "winter" ...

## Resetting variable names to make it easier to reshape dataset after adding moving averages
temp1 <- temp[,5:9]
colnames(temp1) <- paste(colnames(temp1), "raw", sep = ".")
temp <- cbind(temp[,c(1,2,3,4,10,11)], temp1)

Now lets smooth out some the the effect of daily variation and get some moving averages of temporal trends. We will use a centred moving average to capture coldest and hottest periods in each year:

library(tidyverse)
temp11 <- temp %>%
  select(date, date1, yearday, Year, month, season, Delma.raw, RasGhanada.raw, Saadiyat.raw, AlAqah.raw, Dibba.raw) %>%
  mutate(Delma.CMA = rollmean(Delma.raw, k = 61, fill = NA),
         RasGhanada.CMA = rollmean(RasGhanada.raw, k = 61, fill = NA),
         Saadiyat.CMA = rollmean(Saadiyat.raw, k = 61, fill = NA),
         AlAqah.CMA = rollmean(AlAqah.raw, k = 61, fill = NA),
         Dibba.CMA = rollmean(Dibba.raw, k = 61, fill = NA))

## Reformatting and making some plots to visualize SST. Only looking at SST from one site and two different years:
##
temp11$date1 <- as.Date(temp11$date1, '%d/%m/%y')
str(temp)

## 'data.frame':    12370 obs. of  11 variables:
##  $ date          : chr  "1985JAN01" "1985JAN02" "1985JAN03" "1985JAN04" ...
##  $ date1         : chr  "1/1/85" "2/1/85" "3/1/85" "4/1/85" ...
##  $ yearday       : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Year          : int  1985 1985 1985 1985 1985 1985 1985 1985 1985 1985 ...
##  $ month         : chr  "JAN" "JAN" "JAN" "JAN" ...
##  $ season        : chr  "winter" "winter" "winter" "winter" ...
##  $ Delma.raw     : num  22.7 22.6 22.5 22.2 22.1 ...
##  $ RasGhanada.raw: num  23.3 23.1 23 22.8 22.8 ...
##  $ Saadiyat.raw  : num  23.3 23.1 23 22.7 22.7 ...
##  $ AlAqah.raw    : num  24 23.9 23.9 23.9 23.9 ...
##  $ Dibba.raw     : num  24 23.9 23.9 23.8 23.8 ...

require(ggplot2)
temp98 <- subset(temp11, Year ==2002 | Year ==2003)
ggplot( data = temp98, aes(date1, Delma.raw )) + geom_point() + 
  #geom_line(data = temp98, aes(date1, Delma.raw ), color = "red")+
  geom_line(data = temp98, aes(date1, Delma.CMA ), color = "blue")

The data needs reshaping so raw SST (Delma, RasGhanada, Saadiyat, AlAqah, Dibba) and 61 day centered averaged temps are in single columns with a site factor level instead (removing first 30 & last 30 days of data first to remove all NAs produced from the centred moving average)

temp11 <- temp11[31:(nrow(temp11)-31),]
templong <- temp11 %>% 
  gather("key", "value", -date, -date1, -yearday, -Year, -month, -season) %>% 
  separate(key, c("Site", "raw_CMA"), sep = "\\.") %>% 
  spread(raw_CMA, value)

xyplot(raw ~ yearday | as.factor(Year) , data = templong,
       group = Site, t="l",
       grid = TRUE, auto.key=list(space="top", columns=5, 
                                  title="Site", cex.title=1,
                                  lines=TRUE, points=FALSE))

So we can see the expected annually cyclical nature of SST here.

Summarizing the SST data into relevent yearly summaries

Because we are interested in summarizing variables for each calendar year, we need a way to discriminate the different seasons from the temperature data itself in order to summarize, say seasonal(summer / winter) average temperatures or to use the rates of change between summer highs to winter and winter lows to summer each year. We need to do this because using static calendar dates to define these may not capture the appropriate data range for particular summary, because of temporal variation in seasons (e,g, a “late summer” or “early winter”. To begine with, lets detect and plot the90 days after coldest (blue) and hottest day of each year are detected (red) - checking a few different sites & only looking for a few different years 2003-2008

plot_df <- templong %>% 
  mutate(year = lubridate::year(date1)) %>%
  group_by(year, Site) %>%
  mutate(raw_min = +(raw == min(raw)),
         raw_max = +(raw == max(raw))) %>%
  ungroup() %>%
  mutate(raw_min = cumsum(raw_min - lag(raw_min, 90, default = 0)),
         raw_max = cumsum(raw_max - lag(raw_max, 90, default = 0)))
ggplot(
  plot_df %>%
    filter(Site == "Delma", Year >= 2003, Year <= 2008),
  aes(date1, raw)
) +
  geom_line() +
  geom_vline(
    aes(xintercept = date1),
    data = plot_df %>% filter(Site == "Delma", Year >= 2003, Year <= 2008, raw_min > 0),
    alpha = 0.1,
    colour = "blue"
  ) +
  geom_vline(
    aes(xintercept = date1),
    data = plot_df %>% filter(Site == "Delma", Year >= 2003, Year <= 2008, raw_max > 0),
    alpha = 0.1,
    colour = "red"
  )

ggplot(
  plot_df %>%
    filter(Site == "Dibba", Year >= 2003, Year <= 2008),
  aes(date1, raw)
) +
  geom_line() +
  geom_vline(
    aes(xintercept = date1),
    data = plot_df %>% filter(Site == "Dibba", Year >= 2003, Year <= 2008, raw_min > 0),
    alpha = 0.1,
    colour = "blue"
  ) +
  geom_vline(
    aes(xintercept = date1),
    data = plot_df %>% filter(Site == "Dibba", Year >= 2003, Year <= 2008, raw_max > 0),
    alpha = 0.1,
    colour = "red"
  )

We can now subset the raw data to be from the yearly minimum SST at each site for the 90 days after.

raw_minroc <- subset(plot_df, raw_min > 0)
xyplot(raw ~ yearday | as.factor(Year) , data = raw_minroc,
       group = Site, type ="l", xlim = c(0,200),
       grid = TRUE, auto.key=list(space="top", columns=5, 
                                  title="Site", cex.title=1,
                                  lines=TRUE, points=FALSE))

Ideally if we fit a simple linear regression to this line, it would give a measure of the rate of change toward summer from the coldest day of winter in that year. But it seems like a few problems to deal with first - For example, 2005 something weird has gone on at the site RasG. Other instances is where the coldestr day is on the 31st December, (RasG and Saadiyat) which bugs up the by year filtering, so we need to remove (RasG and Saadiyat, 31 dec 2004) running filtering again without those rows

templong <- templong %>%
  filter(date1 != "2004-12-31" | Site != "RasGhanada" ) %>%
  filter(date1 != "2004-12-31" |  Site != "Saadiyat")
nrow(templong)

## [1] 61543

plot_df <- templong %>% 
  mutate(year = lubridate::year(date1)) %>%
  group_by(year, Site) %>%
  mutate(raw_min = +(raw == min(raw)),
         raw_max = +(raw == max(raw))) %>%
  ungroup() %>%
  mutate(raw_min = cumsum(raw_min - lag(raw_min, 90, default = 0)),
         raw_max = cumsum(raw_max - lag(raw_max, 90, default = 0)))
str(plot_df)

## tibble [61,543 × 12] (S3: tbl_df/tbl/data.frame)
##  $ date   : chr [1:61543] "1985JAN31" "1985FEB01" "1985FEB02" "1985FEB03" ...
##  $ date1  : Date[1:61543], format: "1985-01-31" "1985-02-01" ...
##  $ yearday: int [1:61543] 31 32 33 34 35 36 37 38 39 40 ...
##  $ Year   : int [1:61543] 1985 1985 1985 1985 1985 1985 1985 1985 1985 1985 ...
##  $ month  : chr [1:61543] "JAN" "FEB" "FEB" "FEB" ...
##  $ season : chr [1:61543] "winter" "winter" "winter" "winter" ...
##  $ Site   : chr [1:61543] "Delma" "Delma" "Delma" "Delma" ...
##  $ CMA    : num [1:61543] 20.9 20.9 20.8 20.8 20.7 ...
##  $ raw    : num [1:61543] 21.1 20.9 21 21 21.2 ...
##  $ year   : num [1:61543] 1985 1985 1985 1985 1985 ...
##  $ raw_min: int [1:61543] 0 0 0 0 0 0 0 0 0 0 ...
##  $ raw_max: int [1:61543] 0 0 0 0 0 0 0 0 0 0 ...

Now lets plotting again to make sure 90 days after coldest (blue) and hottest day of each year are detected (red) - checking a few different sites

ggplot(plot_df %>% filter(Site =="Delma"), aes(date1, raw)) + 
  geom_line() + 
  geom_vline(aes(xintercept = date1), plot_df %>% filter(Site =="Delma", raw_min > 0),
             alpha = 0.1, colour = "blue") +
  geom_vline(aes(xintercept = date1), plot_df %>% filter(Site =="Delma", raw_max > 0),
             alpha = 0.1, colour = "red")

ggplot(plot_df %>% filter(Site =="Dibba"), aes(date1, raw)) + 
  geom_line() + 
  geom_vline(aes(xintercept = date1), plot_df %>% filter(Site =="Dibba", raw_min > 0),
             alpha = 0.1, colour = "blue") +
  geom_vline(aes(xintercept = date1), plot_df %>% filter(Site =="Dibba", raw_max > 0),
             alpha = 0.1, colour = "red")

Now lets check to see if the data are now amenable to generating yearly rates of change between seasons for the 90 days after the coldest and hottest day of each year.

##checking
raw_minroc <- subset(plot_df, raw_min > 0)
raw_maxroc <- subset(plot_df, raw_max > 0)
## winter to summer
xyplot(raw ~ yearday | as.factor(Year) , data = raw_minroc,
       group = Site, type ="l", xlim = c(0,200),
       grid = TRUE, auto.key=list(space="top", columns=5, 
                                  title="Site", cex.title=1,
                                  lines=TRUE, points=FALSE))

##Summer to winter
xyplot(raw ~ yearday | as.factor(Year) , data = raw_maxroc,
       group = Site, type ="l", #xlim = c(0,200),
       grid = TRUE, auto.key=list(space="top", columns=5, 
                                  title="Site", cex.title=1,
                                  lines=TRUE, points=FALSE))

There appears to be a problem related to when the data starts for sites RasG and Al Aqah for summer to winter in 1985 but this will be discarded anyway (since coral data is only from 2000).

In order to gather a measure of seasonal rates of change, we will fit a linear model within each year for SST~yearday where yearday will be normalized to the coldest (or hottest) day to day 1 before fitting regression so that intercepts are normalised to hottest/coldest day (as applicable)

library(data.table)
str(raw_minroc)

## tibble [15,304 × 12] (S3: tbl_df/tbl/data.frame)
##  $ date   : chr [1:15304] "1985MAR08" "1985MAR09" "1985MAR10" "1985MAR11" ...
##  $ date1  : Date[1:15304], format: "1985-03-08" "1985-03-09" ...
##  $ yearday: int [1:15304] 67 68 69 70 71 72 73 74 75 76 ...
##  $ Year   : int [1:15304] 1985 1985 1985 1985 1985 1985 1985 1985 1985 1985 ...
##  $ month  : chr [1:15304] "MAR" "MAR" "MAR" "MAR" ...
##  $ season : chr [1:15304] "spring" "spring" "spring" "spring" ...
##  $ Site   : chr [1:15304] "Delma" "Delma" "Delma" "Delma" ...
##  $ CMA    : num [1:15304] 20.6 20.6 20.6 20.7 20.7 ...
##  $ raw    : num [1:15304] 18.8 18.8 18.9 19 18.9 ...
##  $ year   : num [1:15304] 1985 1985 1985 1985 1985 ...
##  $ raw_min: int [1:15304] 1 1 1 1 1 1 1 1 1 1 ...
##  $ raw_max: int [1:15304] 0 0 0 0 0 0 0 0 0 0 ...

str(raw_maxroc)

## tibble [15,281 × 12] (S3: tbl_df/tbl/data.frame)
##  $ date   : chr [1:15281] "1985SEP13" "1985SEP14" "1985SEP15" "1985SEP16" ...
##  $ date1  : Date[1:15281], format: "1985-09-13" "1985-09-14" ...
##  $ yearday: int [1:15281] 256 257 258 259 260 261 262 263 264 265 ...
##  $ Year   : int [1:15281] 1985 1985 1985 1985 1985 1985 1985 1985 1985 1985 ...
##  $ month  : chr [1:15281] "SEP" "SEP" "SEP" "SEP" ...
##  $ season : chr [1:15281] "fall" "fall" "fall" "fall" ...
##  $ Site   : chr [1:15281] "Delma" "Delma" "Delma" "Delma" ...
##  $ CMA    : num [1:15281] 33.3 33.3 33.2 33.2 33.2 ...
##  $ raw    : num [1:15281] 35.7 35 34.4 34.6 34.4 ...
##  $ year   : num [1:15281] 1985 1985 1985 1985 1985 ...
##  $ raw_min: int [1:15281] 0 0 0 0 0 0 0 0 0 0 ...
##  $ raw_max: int [1:15281] 1 1 1 1 1 1 1 1 1 1 ...

raw_minroc$Site <- as.factor(raw_minroc$Site)
raw_maxroc$Site <- as.factor(raw_maxroc$Site)

raw_minroc <- raw_minroc %>%
  group_by(Site, year) %>%
  arrange(date1) %>%
  mutate(normday = yearday - min(yearday) + 1)

raw_maxroc <- raw_maxroc %>%
  group_by(Site, year) %>%
  arrange(date1) %>%
  mutate(normday = yearday - min(yearday) + 1)

## 

library(dplyr)
library(tidyr)
library(purrr)
library(broom)

# Fit winter models by Site and year for SST rates of change:

wintrocfin <- raw_minroc %>%
  group_by(Site, year) %>%
  nest() %>%
  mutate(
    model = map(data, ~ lm(raw ~ normday, data = .x)),
    tidied = map(model, tidy),
    glanced = map(model, glance)
  ) %>%
  unnest(cols = c(tidied), names_sep = "_") %>%
  filter(tidied_term == "normday") %>%
  mutate(
    winSlope = tidied_estimate,
    winIntercept = map_dbl(model, ~ coef(.x)[1]),
    winRSquared = map_dbl(glanced, "r.squared")
  ) %>%
  select(Site, year, winRSquared, winIntercept, winSlope) %>%
  left_join(
    raw_minroc %>%
      group_by(Site, year) %>%
      summarise(ccoldday = min(yearday), .groups = "drop"),
    by = c("Site", "year")
  )

# Fit summer models by Site and year for SST rates of change:

summrocfin <- raw_maxroc %>%
  group_by(Site, year) %>%
  nest() %>%
  mutate(
    model = map(data, ~ lm(raw ~ normday, data = .x)),
    tidied = map(model, tidy),
    glanced = map(model, glance)
  ) %>%
  unnest(cols = c(tidied), names_sep = "_") %>%
  filter(tidied_term == "normday") %>%
  mutate(
    sumSlope = tidied_estimate,
    sumIntercept = map_dbl(model, ~ coef(.x)[1]),
    sumRSquared = map_dbl(glanced, "r.squared")
  ) %>%
  select(Site, year, sumRSquared, sumIntercept, sumSlope) %>%
  left_join(
    raw_maxroc %>%
      group_by(Site, year) %>%
      summarise(warmday = max(yearday), .groups = "drop"),
    by = c("Site", "year")
  )

pairs(summrocfin[,3:6])

Now that we have temperature rates of changed summarized, lets move onto to summarizing the hottest and coldest 61 days of each year to indicate the relative “hotness” of the annual summer and “coldness” of the annual winter at each site. For this we will use moving averages rather than daily temperatures, so as to remove daily noise from finding the hottest and coldest periods of each year.

library(dplyr)
library(lubridate)
library(purrr)

# Adjust year so that DECEMBER is counted as part of the following winter - this is for years where the coldest period is close to the start of the year and therefore an attempt to extract the previous 30 days from the coldest moving average value in each yea is not possible. 

templong <- templong %>%
  mutate(yearAD = ifelse(month == "DEC", Year + 1, Year))

# Identify coldest (min CMA) and hottest (max CMA) days per site-year
sumwin_df <- templong %>%
  group_by(yearAD, Site) %>%
  mutate(
    CMA_min = as.integer(CMA == min(CMA, na.rm = TRUE)),
    CMA_max = as.integer(CMA == max(CMA, na.rm = TRUE))
  ) %>%
  ungroup()

# Helper function to extract 61-day window around target indices
extract_window <- function(df, flag_col) {
  inds <- which(df[[flag_col]] == 1)
  rows <- unique(unlist(lapply(inds, function(i) (i - 30):(i + 30))))
  rows <- rows[rows > 0 & rows <= nrow(df)]  # Keep only valid indices
  result <- df[rows, ]
  result[[flag_col]] <- 1  # Mark entire window with flag
  return(result)
}

# Extract 61-day windows for coldest (winter) and hottest (summer) days
sumwin_win <- extract_window(sumwin_df, "CMA_min") %>%
  mutate(yearAD = ifelse(month == "DEC", Year + 1, Year))  # Ensure DEC counted in correct year

sumwin_sum <- extract_window(sumwin_df, "CMA_max")  # Summer doesn't need year fix

# Merge window flags back into full dataset
temp1 <- templong %>%
  left_join(sumwin_sum %>% select(date1, Site, CMA_max), by = c("date1", "Site")) %>%
  left_join(sumwin_win %>% select(date1, Site, CMA_min), by = c("date1", "Site"))

Lets now have a look at a few sites and time period to ensure we are indeed extracting the hottest and coldest periods of each season in each year at each site

We can also visualise varition between sites

Looks good. Now we can easily summarize to extract minimum, maximum, range, mean and the SD (standard deviation) in temperature (representing temperature variability) for the hottest and coldest 61 days of each year at each site.

library(dplyr)
library(purrr)

# Create winteryear column for situations where the coldest period overlaps with the start of the year
temp1 <- temp1 %>%
  mutate(
    Site = as.factor(Site),
    winteryear = ifelse(month == "DEC", Year + 1, Year)
  )

# Subset winter and summer temperature windows
wtemp <- temp1 %>% filter(CMA_min == 1)
stemp <- temp1 %>% filter(CMA_max == 1)

# Generic summarisation function for SST data
summarise_sst <- function(data, group_var, prefix) {
  data %>%
    group_by(across(all_of(group_var)), Site) %>%
    summarise(
      maxM = max(CMA, na.rm = TRUE),
      minM = min(CMA, na.rm = TRUE),
      avM  = mean(CMA, na.rm = TRUE),
      sdM  = sd(CMA, na.rm = TRUE),
      maxR = max(raw, na.rm = TRUE),
      minR = min(raw, na.rm = TRUE),
      avR  = mean(raw, na.rm = TRUE),
      sdR  = sd(raw, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    mutate(site = Site) %>%
    select(-Site) %>%
    rename_with(~ paste0(prefix, .), -c(all_of(group_var), site))
}

# Get summaries
Wtsumall <- summarise_sst(wtemp, "winteryear", "Wt") %>%
  rename(year = winteryear) %>%
  mutate(
    WtrangeM = WtmaxM - WtminM,
    WtrangeR = WtmaxR - WtminR
  ) %>%
  drop_na()

Stsumall <- summarise_sst(stemp, "Year", "St") %>%
  rename(year = Year) %>%
  mutate(
    StrangeM = StmaxM - StminM,
    StrangeR = StmaxR - StminR
  ) %>%
  drop_na()

Finally lets combine these temperature summaries with the coral growth data:

# Prepare coral growth dataset
onyr <- onyr %>%
  mutate(
    site = as.factor(site),
    Site = site
  )

# Rename 'site' columns in SST rate-of-change summaries
wintrocfin <- wintrocfin %>%
  rename(site = Site)

summrocfin <- summrocfin %>%
  rename(site = Site)

# Merge everything
test <- onyr %>%
  left_join(Wtsumall, by = c("year", "site")) %>%
  left_join(Stsumall, by = c("year", "site")) %>%
  left_join(wintrocfin, by = c("year", "site")) %>%
  left_join(summrocfin, by = c("year", "site"))

str(test)

## 'data.frame':    768 obs. of  41 variables:
##  $ year        : num  2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 ...
##  $ site        : Factor w/ 5 levels "AlAqah","Delma",..: 1 1 3 3 3 3 4 4 4 4 ...
##  $ colony      : chr  "ala12" "ala12" "dib01" "dib01" ...
##  $ track       : chr  "t1" "t2" "t1" "t2" ...
##  $ sample      : chr  "ala12_t1" "ala12_t2" "dib01_t1" "dib01_t2" ...
##  $ calc        : num  0.797 1.007 0.956 0.526 0.94 ...
##  $ extension   : num  0.656 0.57 0.747 0.378 0.817 ...
##  $ density     : num  1.53 1.4 1.28 1.39 1.15 ...
##  $ hdbands     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ partmort    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ height      : num  13.5 13.5 NA NA 14 14 7.5 7.5 12 12 ...
##  $ dia         : int  23 23 NA NA 28 28 14 14 16 16 ...
##  $ Site        : Factor w/ 5 levels "AlAqah","Delma",..: 1 1 3 3 3 3 4 4 4 4 ...
##  $ WtmaxM      : num  24.1 24.1 24 24 24 ...
##  $ WtminM      : num  23.1 23.1 23 23 23 ...
##  $ WtavM       : num  23.4 23.4 23.3 23.3 23.3 ...
##  $ WtsdM       : num  0.244 0.244 0.229 0.229 0.229 ...
##  $ WtmaxR      : num  24.1 24.1 24 24 24 ...
##  $ WtminR      : num  22.4 22.4 22.3 22.3 22.3 ...
##  $ WtavR       : num  23.1 23.1 23 23 23 ...
##  $ WtsdR       : num  0.424 0.424 0.417 0.417 0.417 ...
##  $ WtrangeM    : num  0.994 0.994 0.908 0.908 0.908 ...
##  $ WtrangeR    : num  1.76 1.76 1.75 1.75 1.75 ...
##  $ StmaxM      : num  32.3 32.3 32.3 32.3 32.3 ...
##  $ StminM      : num  31.7 31.7 31.8 31.8 31.8 ...
##  $ StavM       : num  32.2 32.2 32.1 32.1 32.1 ...
##  $ StsdM       : num  0.159 0.159 0.137 0.137 0.137 ...
##  $ StmaxR      : num  33.4 33.4 33.3 33.3 33.3 ...
##  $ StminR      : num  31 31 31.6 31.6 31.6 ...
##  $ StavR       : num  32.3 32.3 32.3 32.3 32.3 ...
##  $ StsdR       : num  0.545 0.545 0.432 0.432 0.432 ...
##  $ StrangeM    : num  0.634 0.634 0.534 0.534 0.534 ...
##  $ StrangeR    : num  2.39 2.39 1.71 1.71 1.71 ...
##  $ winRSquared : num  0.874 0.874 0.865 0.865 0.865 ...
##  $ winIntercept: num  21 21 20.9 20.9 20.9 ...
##  $ winSlope    : num  0.0959 0.0959 0.0938 0.0938 0.0938 ...
##  $ ccoldday    : int  42 42 42 42 42 42 42 42 42 42 ...
##  $ sumRSquared : num  0.717 0.717 0.759 0.759 0.759 ...
##  $ sumIntercept: num  32.5 32.5 32.6 32.6 32.6 ...
##  $ sumSlope    : num  -0.0338 -0.0338 -0.0371 -0.0371 -0.0371 ...
##  $ warmday     : int  301 301 301 301 301 301 290 290 290 290 ...

DHW (Degree heating weeks) & DCW (Degree cooling weeks)

Degree heating weeks is a measure of the cumulative effect of temperature (originally heat stress, but optionally cold stress) developed by NOAAs coral reef watch, with details here.

To calculate these metrics we will use the package dbcaDHW. DHW is calculated by comparing temperatures to a climatological threshold - The Maximum Monthly Mean (MMM) for each Site, which is calculated by averaging the annual max monthly means over a baseline climatology period (1985–1995).

# Read and prep data
temp <- read.csv("rawdata/ctemp_sst_mod.csv")
str(temp)

## 'data.frame':    12370 obs. of  11 variables:
##  $ date      : chr  "1985JAN01" "1985JAN02" "1985JAN03" "1985JAN04" ...
##  $ date1     : chr  "1/1/85" "2/1/85" "3/1/85" "4/1/85" ...
##  $ yearday   : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Year      : int  1985 1985 1985 1985 1985 1985 1985 1985 1985 1985 ...
##  $ Delma     : num  22.7 22.6 22.5 22.2 22.1 ...
##  $ RasGhanada: num  23.3 23.1 23 22.8 22.8 ...
##  $ Saadiyat  : num  23.3 23.1 23 22.7 22.7 ...
##  $ AlAqah    : num  24 23.9 23.9 23.9 23.9 ...
##  $ Dibba     : num  24 23.9 23.9 23.8 23.8 ...
##  $ month     : chr  "JAN" "JAN" "JAN" "JAN" ...
##  $ season    : chr  "winter" "winter" "winter" "winter" ...

sites <- names(temp)[5:9]
g_num <- length(sites) + 1
temp1 <- temp[,c(2,5:9)]
temp1$date1 <- as.Date(temp1$date1, '%d/%m/%y')

## Note this is MODIFIED FROM ORIGINAL DHW FORMULA
## create a mean of the maximum monthly mean each year at each site from "past period"
mmtdata <- temp1 %>%
  tidyr::gather("Site", "sst", 2:all_of(g_num)) %>%
  #dplyr::filter(date1 < "2013-12-31") %>% #define past
  dplyr::filter(date1 >= "1985-01-01" & date1 <= "1995-12-31") %>% #define past
  dplyr::mutate(month = month(date1),
                year = year(date1)) %>%
  dplyr::group_by(Site, year, month) %>%
  dplyr::summarise(avg = mean(sst)) %>% #obtain mthly avgs
  dplyr::group_by(Site, year) %>%
  dplyr::summarise(mmmt = max(avg)) %>% #obtain max mthly avg each year over this period
  dplyr::group_by(Site) %>%
  dplyr::summarise(mmmt = mean(mmmt)) #average this value as threshold for DHW

# The DHW calculation is assisted by creating a function to do the work of a rolling 12 week “sum”. 
# A present period also needs to be defined.

## function for rolling cusum of 12 weeks - 84 days)
dhw_calc_12 <- tibbletime::rollify(sum, window = 84)

## observation period, hotspot & dhw calc
dhw <- temp1 %>%
  tidyr::gather("Site", "sst", 2:all_of(g_num)) %>%
  dplyr::filter(date1 >= "2000-01-01" & date1 <= "2013-12-31") %>% #define present
  dplyr::mutate(month = month(date1)) %>%
  dplyr::full_join(mmtdata, by = "Site") %>%
  dplyr::mutate(hspt = sst - mmmt,
                hsptm = ifelse(hspt >= 1, hspt, 0),
                dhw = dhw_calc_12(hsptm)*(1/7)) %>% #hotspot based on greater than 1 degree as per coral watch 
  dplyr::select(-month, -hsptm)

dhwm <- dhw %>%
  dplyr::mutate(year = year(date1)) %>%
  dplyr::group_by(Site, year) %>%
  dplyr::summarise(cummu_dhwm = sum(dhw, na.rm = TRUE))

The dhw object contains daily DWH values for each site. Because we need a summary of dhw by year, we can simply sum all DHW within each year as a yearly metric of total heat stress that year.

dhwm <- dhw %>%
  dplyr::mutate(year = year(date1)) %>%
  dplyr::group_by(Site, year) %>%
  dplyr::summarise(cummu_dhwm = sum(dhw, na.rm = TRUE))

Here is what the DHW data look like for each site over the coral growth period at each site with NOAA bleaching alert levels indicated:

# Converting Site to factor for ordering
dhw <- dhw %>%
  mutate(Site = as.factor(Site))

# Plot: DHW time series with alert thresholds
ggplot(dhw, aes(x = date1, y = dhw)) +
  geom_line(color = "darkorange") +
  geom_hline(yintercept = 4, linetype = "dashed", color = "blue") +
  geom_hline(yintercept = 8, linetype = "dashed", color = "red") +
  facet_wrap(~ Site, scales = "free_y", ncol = 1) +
  labs(
    title = "Degree Heating Weeks (DHW) 2000–2013",
    subtitle = "Dashed lines: Bleaching Alert Levels (4°C-weeks = Blue, 8°C-weeks = Red)",
    x = "Date",
    y = "DHW (°C-weeks)"
  ) +
  theme_minimal()

There is not a lot of variation in dhw, with most daily values throughout most years being zero. This is because for DWH at least, the reference period defining the threshold had generally greater temperatures compared to the period between 2000 and 2013.

The plot below shows the mean monthly temperature (MMT) within each year at each site. The red dashed line indicates the the maximum MMT and the blue dashed line the minimum MMT across all years at that site. The filled lines are the lowess fits illustrating the temporal trend in the data for the maximum MMT (red) and minimum MMT (blue). There are few dhw values > 0 in our dataset because there is a general trend that years (2000-2013) experienced lower maxmimum MMT compared to the reference period (1985-1995).

This plot also illustrates that in winter there is not as clear a temporal signal in minimum MMT, but there is arguably a small downtrend at a couple of sites.

Lets also calculate degree cooling weeks (DCW), where like DHW, the accumulated cold stress over a rolling 12-week (84-day) window is calculated, where daily sea surface temperature (SST) is at least 1°C below a site-specific cold threshold, and the excess cold is summed and expressed in °C-weeks

mmtdata_cool <- temp1 %>%
  tidyr::gather("Site", "sst", 2:all_of(g_num)) %>%
  #dplyr::filter(date1 < "2013-12-31") %>% #define past
  dplyr::filter(date1 >= "1985-01-01" & date1 <= "1995-12-31") %>% #define past
  dplyr::mutate(month = month(date1),
                year = year(date1)) %>%
  dplyr::group_by(Site, year, month) %>%
  dplyr::summarise(avg = mean(sst)) %>% #obtain mthly avgs
  dplyr::group_by(Site, year) %>%
  dplyr::summarise(mmmt = min(avg)) %>% #obtain max mthly avg each year over this period
  dplyr::group_by(Site) %>%
  dplyr::summarise(mmmt = mean(mmmt)) #average this value as threshold for DHW

dhw_calc_12 <- tibbletime::rollify(sum, window = 84)

## observation period, hotspot & dhw calc
dhw_cool <- temp1 %>%
  tidyr::gather("Site", "sst", 2:all_of(g_num)) %>%
  dplyr::filter(date1 >= "2000-01-01" & date1 <= "2013-12-31") %>% #define present
  dplyr::mutate(month = month(date1)) %>%
  dplyr::full_join(mmtdata_cool, by = "Site") %>%
  dplyr::mutate(hspt = sst - mmmt,
                hsptm = ifelse(hspt <= -1, hspt, 0),
                dcw = dhw_calc_12(hsptm)*(1/7)) %>% #hotspot based on less than 1 degree
  dplyr::select(-month, -hsptm)

## again, because we need a summary of dhw by year, we can simply sum all DHW within each year as a yearly metric of total heat stress that year.
dhwm_cool <- dhw_cool %>%
  dplyr::mutate(year = year(date1)) %>%
  dplyr::group_by(Site, year) %>%
  dplyr::summarise(cummu_dcwm = sum(dcw, na.rm = TRUE))

And a plot of DCW values over time: